Morphology of ruff sperm
Martin Bulla
2022-05-02 13:19:56
Code to load tools and data
# TOOLS
require(here)
source(here::here('R/tools.R'))
require(arm)
require(effects)
require(ggpubr)
require(ggsci)
require(grid)
require(gridExtra)
require(magrittr)
require(multcomp)
require(PerformanceAnalytics)
require(rptR)
require(stringi)
require(viridis)
# constants
round_ = 3 # number of decimal places to round model coefficients
nsim = 5000 # number of simulations to extract estimates and 95%CrI
ax_lines = "grey60" # defines color of the axis lines
colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
# functions
getime = function (x) {ifelse(is.na(x), as.numeric(NA), as.numeric(difftime(x, trunc(x,"day"), units = "hours")))}
getDay = function (x) {as.Date(trunc(x, "day"))}
# for R output based on sim
R_out = function(name = "define", model = m, nsim = 5000){
bsim <- sim(model, n.sim=nsim)
l=data.frame(summary(model)$varcor)
l = l[is.na(l$var2),]
l$var1 = ifelse(is.na(l$var1),"",l$var1)
l$pred = paste(l$grp,l$var1)
q50={}
q025={}
q975={}
pred={}
# variance of random effects
for (ran in names(bsim@ranef)) {
#ran =names(bsim@ranef)[1]
ran_type = l$var1[l$grp == ran]
for(i in ran_type){
# i = ran_type[2]
q50=c(q50,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.5)))
q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.025)))
q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.975)))
pred= c(pred,paste(ran, i))
}
}
# residual variance
q50=c(q50,quantile(bsim@sigma^2, prob=c(0.5)))
q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
pred= c(pred,'Residual')
ci = c(round(100*q025/sum(q025))[1], round(100*q975/sum(q975))[1])
ci = ci[order(ci)]
ri=data.table(model = name, repeatability=paste0(round(100*q50/sum(q50)),'%')[1], CI = paste0(paste(ci[1], ci[2], sep ="-"), '%'))
return(ri)
}
# DATA
# composite measures
x = fread(here::here('Data/DAT_morpho.csv'))
setnames(x,old = 'pic', new = 'sperm_ID')
x[, sample_ID:=as.character(sample_ID)]
# mneasurements from minip vs rest are the same - so no need to control for
#ggplot(x[part == 'Tail'], aes(x = manip, y = Pixels)) + geom_boxplot()
d1 = x[,.(bird_ID, sample_ID, sperm_ID, part, Pixels)]
d2 = d1[, sum(Pixels), by = list(bird_ID, sample_ID, sperm_ID)]
names(d2)[4] = 'Pixels'
d2[,part := 'Total']
d3 = d1[part%in%c('Midpiece', 'Tail'), sum(Pixels), by = list(bird_ID, sample_ID, sperm_ID)]
names(d3)[4] = 'Pixels'
d3[,part := 'Flagellum']
d4 = d1[part%in%c('Acrosome', 'Nucleus'), sum(Pixels), by = list(bird_ID, sample_ID, sperm_ID)]
names(d4)[4] = 'Pixels'
d4[,part := 'Head']
b = merge(d1,d2, all = TRUE)
b = merge(b,d3, all = TRUE)
b = merge(b,d4, all = TRUE)
b[, Length_µm:=Pixels*0.078]
# add metadata and motility
b = merge(b, x[,.(bird_ID, sample_ID, sperm_ID, part, manip)], all.x = TRUE, by=c('bird_ID', 'sample_ID', 'sperm_ID', 'part'))
b[manip%in%"", manip:=NA]
d = data.table(read_excel(here::here('Data/motility.xlsx'), sheet = 1))
s = data.table(read_excel(here::here('Data/sampling_2021_cleaned.xlsx'), sheet = 1))
s = s[!is.na(recording)]
m = data.table(read_excel(here::here('Data/ruff_males_Seewiesen.xlsx'), sheet = 1))#, range = "A1:G161"))
m[, hatch_year:=as.numeric(substr(Hatch_date,1,4)) ]
m[, age := 2021-hatch_year]
m = m[Ind_ID %in%unique(s$DB_ID[!is.na(s$DB_ID)]),.(Ind_ID, Morph, age)]
s = merge(s,m, by.x = 'DB_ID', by.y = 'Ind_ID', all.x = TRUE)
b = merge(b, s[,.(sample_ID, bird_ID, Morph, age, datetime, type, sperm, recording, rec_measured, month)], by.x = c('sample_ID', 'bird_ID'), by.y = c('sample_ID', 'bird_ID'), all.x = TRUE)
b = merge(b, d, by.x = c('bird_ID', 'month'), by.y = c('ID', 'date'), all.x = TRUE)
#d = (d[!duplicated(d)]) # shouldn't be necessary, but for some reason the above call duplicates the dataset
b[is.na(Morph), Morph := 'Zebra finch']
b[Morph == 'F', Morph := 'Faeder']
b[Morph == 'I', Morph := 'Independent']
b[Morph == 'S', Morph := 'Satellite']
b[is.na(issues), issues := 'zero']
b = b[!Morph %in% 'Zebra finch']
b[, part := factor(part, levels = c('Acrosome', 'Nucleus', 'Midpiece', 'Tail','Head','Flagellum', 'Total'))]
b[part %in% c('Acrosome', 'Nucleus', 'Midpiece', 'Tail'), measure := 'original']
b[part %in% c('Head','Flagellum', 'Total'), measure := 'composite']
#nrow(d) # N 139
# prepare for correlations and repeatability
bw = reshape(b[,.(bird_ID,Morph, age, datetime, month, sample_ID, sperm_ID, VAP,VSL,VCL, motileCount, part, Length_µm)], idvar = c('bird_ID','Morph','age','datetime', 'month', 'sample_ID', 'sperm_ID','VAP','VSL','VCL', 'motileCount'), timevar = 'part', direction = "wide")
setnames(bw,old = c('Length_µm.Acrosome', 'Length_µm.Nucleus','Length_µm.Flagellum','Length_µm.Head','Length_µm.Midpiece','Length_µm.Tail', 'Length_µm.Total'), new = c('Acrosome', 'Flagellum', 'Head','Nucleus', 'Midpiece', 'Tail','Total'))
# add relative measures
bw[, Midpiece_rel := Midpiece/Total]
bw[, Flagellum_rel := Flagellum/Total]
# mean/male dataset
a = b[, list(mean(Length_µm), mean(VAP),mean(VSL), mean(VCL), mean(motileCount)), by = list(bird_ID, Morph, age, part)]
setnames(a, old = c('V1', 'V2','V3','V4','V5'), new = c('Length_avg', 'VAP', 'VSL', 'VCL', 'motileCount'))
a1 = data.table(bird_ID = unique(a$bird_ID), blind = c(rep(c('Independent','Satellite', 'Faeder'), floor(length(unique(a$bird_ID))/3)), 'Independent','Satellite'))
a = merge(a,a1, all.x = TRUE)
aw = reshape(a, idvar = c('bird_ID','Morph', 'blind','age','VAP','VSL','VCL', 'motileCount'), timevar = "part", direction = "wide")
names(aw) = c('bird_ID','Morph', 'VAP','VSL','VCL', 'motileCount','blind','age',as.character(unique(a$part)))
# add relative measures
aw[, Midpiece_rel := Midpiece/Total]
aw[, Flagellum_rel := Flagellum/Total]
# morph as factor
b[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))]
bw[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))]
a[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))]
aw[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))]
# dataset for relative measurements
bt = b[part == 'Total',.(Morph, bird_ID, sample_ID, sperm_ID, part, Length_µm, VAP, VSL, VCL, motileCount)]
b1 = b[part%in%c('Midpiece'),.(Morph, bird_ID, sample_ID, sperm_ID, part, Length_µm, VAP, VSL, VCL, motileCount)]
b1[, Length_rel := Length_µm/bt$Length_µm]
b2 = b[part%in%c('Flagellum'),.(Morph, bird_ID, sample_ID, sperm_ID, part, Length_µm, VAP, VSL, VCL, motileCount)]
b2[, Length_rel := Length_µm/bt$Length_µm]
br = rbind(b1,b2)
br$Length_µm = NULL
at = a[part == 'Total',.(Morph, bird_ID, part, Length_avg, VAP, VSL, VCL, motileCount)]
a1 = a[part%in%c('Midpiece'),.(Morph, bird_ID, part, Length_avg, VAP, VSL, VCL, motileCount)]
a1[, Length_rel := Length_avg/at$Length_avg]
a2 = a[part%in%c('Flagellum'),.(Morph, bird_ID, part, Length_avg, VAP, VSL, VCL, motileCount)]
a2[, Length_rel := Length_avg/at$Length_avg]
ar = rbind(a1,a2)
ar$Length_avg = NULL
# dataset for correlations
h = b[part == 'Head']
setnames(h, old = "Length_µm", new="Head_µm")
h[, Acrosome_µm := b[part == 'Acrosome',.(Length_µm)]]
h[, Nucleus_µm := b[part == 'Nucleus',.(Length_µm)]]
h[, Midpiece_µm := b[part == 'Midpiece',.(Length_µm)]]
h[, Tail_µm := b[part == 'Tail',.(Length_µm)]]
h[, Flagellum_µm := b[part == 'Flagellum',.(Length_µm)]]
h[, Total_µm := b[part == 'Total',.(Length_µm)]]
h[, Midpiece_rel := br[part == 'Midpiece',.(Length_rel)]]
h[, Flagellum_rel := br[part == 'Flagellum',.(Length_rel)]]
ha = a[part == 'Head']
setnames(ha, old = "Length_avg", new="Head_µm")
ha[, Acrosome_µm := a[part == 'Acrosome',.(Length_avg)]]
ha[, Nucleus_µm := a[part == 'Nucleus',.(Length_avg)]]
ha[, Midpiece_µm := a[part == 'Midpiece',.(Length_avg)]]
ha[, Tail_µm := a[part == 'Tail',.(Length_avg)]]
ha[, Flagellum_µm := a[part == 'Flagellum',.(Length_avg)]]
ha[, Total_µm := a[part == 'Total',.(Length_avg)]]
ha[, Midpiece_rel := ar[part == 'Midpiece',.(Length_rel)]]
ha[, Flagellum_rel := ar[part == 'Flagellum',.(Length_rel)]]
# CV dataset
cv_ = b[, cv(Length_µm), by = list(bird_ID, part, Morph)]
cv_[ , Morph123 := as.numeric(Morph)]
names(cv_) [4]='CV'Exploration
b[, order_ := mean(Length_µm), by = bird_ID]
b_ = b[part =='Total']
b_[,bird_ID := reorder(bird_ID, Length_µm, mean)]
b[, bird_ID := factor(bird_ID, levels = levels(b_$bird_ID))]
g = ggplot(b, aes(x = reorder(as.factor(bird_ID),Length_µm,mean), y = Length_µm)) + facet_wrap(~part, nrow = 7, scales = 'free') + #x = reorder(as.factor(bird_ID),Length_µm,mean)
geom_boxplot(aes(col = Morph)) +
labs(subtitle = 'Distribution within & across males (ordered by total length)')+
#geom_dotplot(binaxis = 'y', stackdir = 'center',
# position = position_dodge(), col = 'red', fill ='red')+
scale_colour_manual(values = colors) +
#scale_color_viridis(discrete=TRUE)+
xlab('Male ID') +
theme_bw() +
theme(axis.text.x = element_blank())
#$, legend.position = "none")
g
ggsave(here::here('Output/morpho_within_male_boxplots_ordered.png'), g, width = 20, height = 15, units = 'cm') #fig.width=5, fig.height = 5
chart.Correlation(bw[, c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
mtext("Single sperm", side=3, line=3)
#dev.copy(png,'Output/VD_corr_single.png')
#dev.off()
# cor_parts_avg, fig.width=5, fig.height=5
chart.Correlation(aw[, c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total','Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
mtext("Male averages", side=3, line=3)
#dev.copy(png,'Output/VD_corr_avg.png')
#dev.off()
!!! On a single sperm level, head strongly correlates with tail and hence both head and tail also strongly correlate with total sperm length. In addition, nucleus strongly correlates with flagellum. !!!
!!! Correlations on ♂ averages differ from those on single sperm, perhaps because some sperm parts do not have a high repeatability (see Repeatability section). Nucleus determines the head length and tail determines flagellum length and also highly correlates with total sperm length. Head (or nucleus) and to a less extreme flagellum (tail) strongly correlate with total sperm length and hence also with relative flagellum length. !!!
Thus, were possible, analyses should be on the level of individual sperm. Whether the correlations are consistent across morphs reveals the next graph.
#fig.width=5, fig.height = 5
chart.Correlation(bw[Morph == 'Independent', c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
mtext("Independent", side=3, line=3)
#dev.copy(png,'Output/morpho_corr_single_I.png')
#dev.off()
#cor_parts_S, fig.width=5, fig.height = 5
chart.Correlation(bw[Morph == 'Satellite', c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
mtext("Satellite", side=3, line=3)
#dev.copy(png,'Output/morpho_corr_single_S.png')
#dev.off()
#cor_parts_F, fig.width=5, fig.height = 5
chart.Correlation(bw[Morph == 'Faeder', c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
mtext("Feader", side=3, line=3)
#dev.copy(png,'Output/morpho_corr_single_F.png')
#dev.off()
!!! Correlations look similar across morphs, except for midpiece in faeders where it more strongly correlates with head, tail and hence also total sperm length !!!
Repeatability
# male
lfsim = list()
lfrpt = list()
for(i in c('Acrosome','Flagellum','Head','Midpiece','Nucleus','Tail','Total')){
part_ = i
# part_ = "Acrosome"
dd = b[part == part_]
m = lmer(Length_µm ~ 1+(1|bird_ID), dd)
Rf = R_out(part_)
lfsim[[i]] = Rf[, method_CI:='arm package']
R = rpt(Length_µm ~ (1 | bird_ID), grname = "bird_ID", data = dd, datatype = "Gaussian")#, nboot = 0, npermut = 0)
RR = data.table(merge(data.frame(name =part_), paste0(round(R$R*100),'%'))) %>% setnames(new = c('part', 'Repeatability'))
RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')]
lfrpt[[i]] = RR[, method_CI := 'rpt package']
print(i)
}
x = do.call(rbind,lfsim)
names(x)[1] = "part"
y = do.call(rbind,lfrpt)
x[, pred:= as.numeric(substr(repeatability,1,2))]
x[, lwr:= as.numeric(substr(CI,1,2))]
x[, upr:= as.numeric(substr(CI,4,5))]
y[, pred:= as.numeric(substr(Repeatability,1,2))]
y[nchar(CI) == 5, CI := paste0(0,CI) ]
y[, lwr:= as.numeric(substr(CI,1,2))]
y[, upr:= as.numeric(substr(CI,4,5))]
names(y)[2] = tolower( names(y)[2])
xy = rbind(x,y)
xy[, part := factor(part, levels=c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total"))]
# morph
lfsim = list()
lfrpt = list()
for(i in c('Acrosome','Flagellum','Head','Midpiece','Nucleus','Tail','Total')){
part_ = i
# part_ = "Acrosome"
dd = b[part == part_]
m = lmer(Length_µm ~ 1+(1|Morph), dd)
Rf = R_out(part_)
lfsim[[i]] = Rf[, method_CI:='arm package']
R = rpt(Length_µm ~ (1 | Morph), grname = "Morph", data = dd, datatype = "Gaussian")#, nboot = 0, npermut = 0)
RR = data.table(merge(data.frame(name =part_), paste0(round(R$R*100),'%'))) %>% setnames(new = c('part', 'Repeatability'))
RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')]
lfrpt[[i]] = RR[, method_CI := 'rpt package']
print(i)
}
x = do.call(rbind,lfsim)
names(x)[1] = "part"
y = do.call(rbind,lfrpt)
x[, pred:= as.numeric(sub("\\%.*", "", repeatability))]
x[, lwr:= as.numeric(sub("\\-.*", "", CI))]
x[, upr:= as.numeric(sub("\\%.*", "",sub(".*-", "", CI)))]
y[, pred:= as.numeric(sub("\\%.*", "", Repeatability))]
y[, lwr:= as.numeric(sub("\\-.*", "", CI))]
y[, upr:= as.numeric(sub("\\%.*", "",sub(".*-", "", CI)))]
#y[nchar(CI) == 5, CI := paste0(0,CI) ]
#y[, lwr:= as.numeric(substr(CI,1,2))]
#y[, upr:= as.numeric(substr(CI,4,5))]
names(y)[2] = tolower( names(y)[2])
xy2 = rbind(x,y)
xy2[, part := factor(part, levels=c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total"))] g1 =
ggplot(xy, aes(x = part, y = pred, col = method_CI)) +
geom_errorbar(aes(ymin = lwr, ymax = upr, col = method_CI), width = 0.1, position = position_dodge(width = 0.25) ) +
#ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.25)) +
scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
labs(x = NULL, y = "Repeatability [%]", subtitle = "within male")+
ylim(c(0,100))+
coord_flip()+
theme_bw() +
theme(plot.title = element_text(size=9),
legend.position = "none")
#ggsave(here::here('Output/morpho_Repeatability_within-male.png'),g1, width = 10, height =7, units = 'cm')
g2 =
ggplot(xy2, aes(x = part, y = pred, col = method_CI)) +
geom_errorbar(aes(ymin = lwr, ymax = upr, col = method_CI), width = 0.1, position = position_dodge(width = 0.25) ) +
#ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.25)) +
scale_color_viridis(discrete=TRUE, begin=0, end = 0.5, guide = guide_legend(reverse = TRUE)) +
#scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
labs(x = NULL, y = "Repeatability [%]", subtitle = "within morph")+
ylim(c(0,100))+
coord_flip()+
theme_bw() +
theme(plot.title = element_text(size=9),
axis.title.y = element_blank(),
axis.text.y = element_blank())
#ggsave(here::here('Output/morpho_Repeatability_within-morph.png'),g, width = 10, height =7, units = 'cm')
grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), size = "last"))
#grid.arrange(g1,g2)
ggsave(here::here('Output/morpho_Repeatability.png'),cbind(ggplotGrob(g1),ggplotGrob(g2) , size = "last"), width = 6*2, height =3*1.5, units = 'cm') Differences - “raw”
Red indicates mean
g1 = # dummy to extract variables for median calculation
ggplot(b, aes(x = Morph, y = Length_µm)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 0.5)+
geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) +
stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
scale_color_viridis(discrete=TRUE, alpha = 0.8)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
ggtitle('Single measurements') +
ylab('Length [µm]') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #,
g2 =
ggplot(a, aes(x = Morph, y = Length_avg)) +
geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 1)+
# position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
scale_color_viridis(discrete=TRUE, alpha = 0.8)+
stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
facet_wrap(~part, scales = 'free_y', nrow = 2)+
ggtitle('Male means') +
ylab('Length [µm]') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #axis.title.x = element_blank(), axis.text.x = element_blank(),
cv_ = b[, cv(Length_µm), by = list(bird_ID, part, Morph)]
cv_[ , Morph123 := as.numeric(Morph)]
names(cv_) [4]='CV'
g3 =
ggplot(cv_, aes(x =Morph , y = CV)) +
geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 1)+
stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
scale_color_viridis(discrete=TRUE, alpha = 0.8)+
facet_wrap(~part, nrow = 2, scales = "free_y") +
guides(x = guide_axis(angle = -45)) +
ylab('Coefficient of variation') +
ggtitle('Male values') +
theme_bw() +
theme(legend.position = "none",
plot.title = element_text(size=9),
axis.title.x = element_blank()
)
gg1 <- ggplotGrob(g1)
gg2 <- ggplotGrob(g2)
gg3 <- ggplotGrob(g3)
grid.draw(rbind(gg1, gg2, gg3, size = "last"))
#grid.arrange(g1,g2)
ggsave(here::here('Output/morpho_boxplots_.png'), rbind(gg1,gg2,gg3, size = "last"), width = 7*3, height =10*1.5, units = 'cm')
#ggsave(here::here('Output/morpho_per_male.png'),g, width = 10, height =10, units = 'cm') g=
ggplot(cv_, aes(x =Morph , y = CV, col = part, fill = part)) +
#geom_dotplot(binaxis = 'y', stackdir = 'center', position = position_dodge(), col = 'grey', dotsize = 0.5)+
geom_boxplot(fill = NA) +
theme_bw() +
theme(axis.title.x = element_blank())
g
ggsave(here::here('Output/CV_per_male_alternative.png'),g, width = 10, height =7, units = 'cm') g1 = # dummy to extract variables for median calculation
ggplot(br, aes(x = Morph, y = Length_rel)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 0.5)+
geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) +
stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
scale_color_viridis(discrete=TRUE, alpha = 0.8)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
ggtitle('Single measurements') +
ylab('Relative to total sperm length') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
#axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #,
g2 = # dummy to extract variables for median calculation
ggplot(ar, aes(x = Morph, y = Length_rel)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 0.5)+
geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) +
stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
scale_color_viridis(discrete=TRUE, alpha = 0.8)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
ggtitle('Male means') +
ylab('Relative to total sperm length') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
plot.title = element_text(size=9)
)
gg1 <- ggplotGrob(g1)
gg2 <- ggplotGrob(g2)
grid.draw(cbind(gg1, gg2, size = "last"))
#grid.arrange(g1,g2)
ggsave(here::here('Output/morpho_boxplots_rel.png'),cbind(gg1,gg2, size = "last"), width = 7*2, height =10*1.5, units = 'cm') Morpho-space
# single sperm
bs = b[, quantile(Length_µm, prob = c(0.5)), by = list(part,Morph)]
names(bs)[3] = 'median'
bs$lwr = b[, quantile(Length_µm, prob = c(0.025)), by = list(part,Morph)]$V1
bs$upr = b[, quantile(Length_µm, prob = c(0.975)), by = list(part,Morph)]$V1
# mean +/-sd
bs = b[, mean(Length_µm), by = list(part,Morph)]
names(bs)[3] = 'median'
bs$lwr = b[, mean(Length_µm)-sd(Length_µm), by = list(part,Morph)]$V1
bs$upr = b[, mean(Length_µm)+sd(Length_µm), by = list(part,Morph)]$V1
bs_f = bs[part == 'Flagellum']
bs_f[,Midpiece_µm:= bs[part == 'Midpiece',.(median)]]
bs_fm = bs[part == 'Midpiece']
bs_fm[,Flagellum_µm:= bs[part == 'Flagellum',.(median)]]
g0 =
ggplot(h, aes(x = Flagellum_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = bs_f, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
geom_segment(data = bs_f, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
geom_segment(data = bs_fm, aes(x = Flagellum_µm, xend = Flagellum_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = bs_f, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#geom_point(data = bs_f, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
#ylim(c(20,28)) +
#xlim(c(24,36)) +
annotate(geom="text", x=80, y=29.5, label="Independent", size = 2, col = colors[1], hjust = 0) +
annotate(geom="text", x=80, y=29, label="Satellite", size = 2, col = colors[2], hjust = 0) +
annotate(geom="text", x=80, y=28.5, label="Faeder", size = 2, col = colors[3], hjust = 0) +
annotate(geom="text", x=80, y=28, label="Mean +/-SD", size = 2, hjust = 0) +
#scale_colour_viridis(discrete=TRUE)+
ggtitle('Single measurements') +
theme_bw() +
theme(
#legend.position=c(.9,.9),
legend.position = "none",
#axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
)
bs_h = bs[part == 'Head']
bs_h[,Midpiece_µm:= bs[part == 'Midpiece',.(median)]]
bs_m = bs[part == 'Midpiece']
bs_m[,Head_µm:= bs[part == 'Head',.(median)]]
g1 =
ggplot(h, aes(x = Head_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = bs_h, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
geom_segment(data = bs_h, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
geom_segment(data = bs_m, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = bs_h, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
#ylim(c(20,28)) +
#xlim(c(24,36)) +
#annotate(geom="text", x=25, y=29.5, label="Independent", size = 2, col = colors[1], hjust = 0) +
#annotate(geom="text", x=25, y=29, label="Satellite", size = 2, col = colors[2], hjust = 0) +
#annotate(geom="text", x=25, y=28.5, label="Faeder", size = 2, col = colors[3], hjust = 0) +
#annotate(geom="text", x=25, y=28, label="Mean +/-SD", size = 2, hjust = 0) +
#scale_colour_viridis(discrete=TRUE)+
ggtitle('Single measurements') +
theme_bw() +
theme(
#legend.position=c(.9,.9),
legend.position = "none",
#axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
)
bs_h_ = bs[part == 'Head']
bs_h_[,Flagellum_µm:= bs[part == 'Flagellum',.(median)]]
bs_f = bs[part == 'Flagellum']
bs_f[,Head_µm:= bs[part == 'Head',.(median)]]
g2 =
ggplot(h, aes(x = Head_µm, y = Flagellum_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(pch = 21)+
geom_point(data = bs_h_, aes(x = median, y = Flagellum_µm), col = 'black', size = 2) +
geom_segment(data = bs_h_, aes(x = lwr, xend = upr, y =Flagellum_µm, yend =Flagellum_µm), col = "black" ) +
geom_segment(data = bs_f, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = bs_h_, aes(x = median, y = Flagellum_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#geom_point(data = bs_h_, aes(x = median, y = Flagellum_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
#scale_colour_viridis(discrete=TRUE)+
theme_bw() +
theme(legend.position = "none",
#axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
)
bs_n = bs[part == 'Nucleus']
bs_n[,Acrosome_µm:= bs[part == 'Acrosome',.(median)]]
bs_a = bs[part == 'Acrosome']
bs_a[,Nucleus_µm:= bs[part == 'Nucleus',.(median)]]
g3 =
ggplot(h, aes(x = Nucleus_µm, y = Acrosome_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = bs_n, aes(x = median, y = Acrosome_µm), col = 'black', size = 2) +
geom_segment(data = bs_n, aes(x = lwr, xend = upr, y =Acrosome_µm, yend =Acrosome_µm), col = "black" ) +
geom_segment(data = bs_a, aes(x = Nucleus_µm, xend = Nucleus_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = bs_n, aes(x = median, y = Acrosome_µm, col = Morph), size = 1) +
# geom_point(data = bs_n, aes(x = median, y = Acrosome_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#scale_colour_viridis(discrete=TRUE)+
theme_bw() +
theme(legend.position = "none",
#axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
)
# male means
as = a[, quantile(Length_avg, prob = c(0.5)), by = list(part,Morph)]
names(as)[3] = 'median'
as$lwr = a[, quantile(Length_avg, prob = c(0.025)), by = list(part,Morph)]$V1
as$upr = a[, quantile(Length_avg, prob = c(0.975)), by = list(part,Morph)]$V1
# mean +/-sd
as = a[, mean(Length_avg), by = list(part,Morph)]
names(as)[3] = 'median'
as$lwr = a[, mean(Length_avg)-sd(Length_avg), by = list(part,Morph)]$V1
as$upr = a[, mean(Length_avg)+sd(Length_avg), by = list(part,Morph)]$V1
as_f = as[part == 'Flagellum']
as_f[,Midpiece_µm:= as[part == 'Midpiece',.(median)]]
as_fm = as[part == 'Midpiece']
as_fm[,Flagellum_µm:= as[part == 'Flagellum',.(median)]]
g0a =
ggplot(ha, aes(x = Flagellum_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = as_f, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
geom_segment(data = as_f, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
geom_segment(data = as_fm, aes(x = Flagellum_µm, xend = Flagellum_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = as_f, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#geom_point(data = aas_f, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
#ylim(c(20,28)) +
#xlim(c(24,36)) +
#scale_colour_viridis(discrete=TRUE)+
ggtitle('Male means') +
theme_bw() +
theme(
#legend.position=c(.9,.9),
legend.position = "none",
#axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
)
as_h = as[part == 'Head']
as_h[,Midpiece_µm:= as[part == 'Midpiece',.(median)]]
as_m = as[part == 'Midpiece']
as_m[,Head_µm:= as[part == 'Head',.(median)]]
g1a =
ggplot(ha, aes(x = Head_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = as_h, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
geom_segment(data = as_h, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
geom_segment(data = as_m, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = as_h, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#geom_point(data = as_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
#ggtitle('Male means') +
theme_bw() +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
plot.title = element_text(size=9)
)
as_h_ = as[part == 'Head']
as_h_[,Flagellum_µm:= as[part == 'Flagellum',.(median)]]
as_f = as[part == 'Flagellum']
as_f[,Head_µm:= as[part == 'Head',.(median)]]
g2a =
ggplot(ha, aes(x = Head_µm, y = Flagellum_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = as_h_, aes(x = median, y = Flagellum_µm), col = 'black', size = 2) +
geom_segment(data = as_h_, aes(x = lwr, xend = upr, y =Flagellum_µm, yend =Flagellum_µm), col = "black" ) +
geom_segment(data = as_f, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = as_h_, aes(x = median, y = Flagellum_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#scale_colour_viridis(discrete=TRUE)+
theme_bw() +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
plot.title = element_text(size=9)
)
as_n = as[part == 'Nucleus']
as_n[,Acrosome_µm:= as[part == 'Acrosome',.(median)]]
as_a = as[part == 'Acrosome']
as_a[,Nucleus_µm:= as[part == 'Nucleus',.(median)]]
g3a =
ggplot(ha, aes(x = Nucleus_µm, y = Acrosome_µm, col = Morph, fill = Morph)) +
geom_point(pch = 21)+
geom_point(data = as_n, aes(x = median, y = Acrosome_µm), col = 'black', size = 2) +
geom_segment(data = as_n, aes(x = lwr, xend = upr, y =Acrosome_µm, yend =Acrosome_µm), col = "black" ) +
geom_segment(data = as_a, aes(x = Nucleus_µm, xend = Nucleus_µm, y =lwr, yend =upr), col = "black" )+
geom_point(data = as_n, aes(x = median, y = Acrosome_µm, col = Morph), size = 1) +
scale_colour_manual(values = colors)+
scale_fill_manual(values = alpha(colors, 0.4))+
#geom_point(data = as_n, aes(x = median, y = Acrosome_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
theme_bw() +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
plot.title = element_text(size=9)
)
# export
gg0 <- ggplotGrob(g0)
gg1 <- ggplotGrob(g1)
gg2 <- ggplotGrob(g2)
gg3 <- ggplotGrob(g3)
gg0a <- ggplotGrob(g0a)
gg1a <- ggplotGrob(g1a)
gg2a <- ggplotGrob(g2a)
gg3a <- ggplotGrob(g3a)
grid.draw(cbind(rbind(gg0,gg1,gg2, gg3, size = "first"), rbind(gg0a,gg1a,gg2a, gg3a, size = "first"), size = "first"))
#grid.arrange(g1,g2)
ggsave(here::here('Output/morpho_corWithMeans.png'),cbind(rbind(gg0,gg1,gg2, gg3, size = "first"), rbind(gg0a,gg1a,gg2a, gg3a, size = "first"), size = "first"), width = 10, height =20, units = 'cm') bw[,morph123 :=ifelse(Morph == 'Independent', 1, ifelse(Morph == 'Satellite', 2,3))]
p = cloud(Tail ~ Head * Midpiece, data=bw, group = Morph, col = colors,
key = list(text = list(c('Independent', 'Satellite','Faeder'),col=colors,cex=0.6)),
pch = 16, cex = 0.8, alpha = 0.75,
xlab=list(cex=0.7),
ylab=list(cex=0.7),
zlab = list(cex=0.7)
)
npanel <- c(4, 2)
rotx <- c(-50, -80)
rotz <- seq(30, 300, length = npanel[1]+1)
update(p[rep(1, prod(npanel))], layout = npanel,
panel = function(..., screen) {
panel.cloud(..., screen = list(z = rotz[current.column()],
x = rotx[current.row()]))
})
png(here::here('Output/morpho_cor_3D.png'),width = 20, height = 10, units = "cm", res = 600)
update(p[rep(1, prod(npanel))], layout = npanel,
panel = function(..., screen) {
panel.cloud(..., screen = list(z = rotz[current.column()],
x = rotx[current.row()]))
})
dev.off()Models
Effect sizes
# morpho
# for averages
l = list()
lp =list()
lpr = list()
lcv = list()
lpcv = list()
for(i in unique(a$part)){
#i ='Nucleus'
m = lm(scale(Length_avg) ~ Morph, a[part == i])
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
l[[i]]=data.frame(response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lm(Length_avg ~ Morph, a[part == i])
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$Length_avg <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$part=i
lp[[i]] = data.table(newD)
print(i)
}
for(i in unique(ar$part)){
if(i == 'Midpiece'){ii = 'Midpiece_rel'}
if(i == 'Flagellum'){ii = 'Flagellum_rel'}
m = lm(scale(Length_rel) ~ Morph, ar[part == i])
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
l[[ii]]=data.frame(response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lm(Length_rel ~ Morph, ar[part == i])
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$Length_rel <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$part=i
lpr[[i]] = data.table(newD)
print(i)
}
for(i in unique(cv_$part)){
#i ='Nucleus'
m = lm(scale(CV) ~ Morph, cv_[part == i])
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
lcv[[i]]=data.frame(response=paste('CV',i),effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lm(CV ~ Morph, cv_[part == i])
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$Length_avg <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$part=i
lpcv[[i]] = data.table(newD)
print(i)
}
ll = data.table(do.call(rbind,l) )
ll[, response := factor(response, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))]
ll[, effect := factor(effect, levels=rev(c("(Intercept)", "MorphSatellite", "MorphFaeder")))]
ll[, unit := 'male average']
llcv = data.table(do.call(rbind,lcv) )
llcv[, response := factor(response, levels=rev(c("CV Acrosome", "CV Nucleus", "CV Head", "CV Midpiece","CV Tail","CV Flagellum","CV Total")))]
llcv[, effect := factor(effect, levels=rev(c("(Intercept)", "MorphSatellite", "MorphFaeder")))]
llp = data.table(do.call(rbind,lp) )
llp[, part := factor(part, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total")))]
llp[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
llpr = data.table(do.call(rbind,lpr) )
llpr[, part := factor(part, levels=rev(c("Midpiece","Flagellum")))]
llpr[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
# for single values
ls = list()
lps = list()
lpsr = list()
for(i in unique(a$part)){
#i ='Nucleus'
m = lmer(scale(Length_µm) ~ Morph + (1|bird_ID), b[part == i])
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=5000)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975))
ls[[i]]=data.frame(response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lmer(Length_µm ~ Morph + (1|bird_ID), b[part == i])
bsim = sim(m, n.sim=nsim)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$Length_µm <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$part=i
lps[[i]] = data.table(newD)
print(i)
}
for(i in unique(ar$part)){
if(i == 'Midpiece'){ii = 'Midpiece_rel'}
if(i == 'Flagellum'){ii = 'Flagellum_rel'}
m = lmer(scale(Length_rel) ~ Morph + (1|bird_ID), br[part == i])
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=5000)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975))
ls[[ii]]=data.frame(response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lmer(Length_rel ~ Morph + (1|bird_ID), br[part == i])
bsim = sim(m, n.sim=nsim)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(ar$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$Length_rel <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$part=i
lpsr[[i]] = data.table(newD)
print(i)
}
lls = data.table(do.call(rbind,ls) )
lls[, response := factor(response, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))]
lls[, effect := factor(effect, levels=rev(c("(Intercept)", "MorphSatellite", "MorphFaeder")))]
lls[, unit := 'single sperm']
llps = data.table(do.call(rbind,lps) )
llps[, part := factor(part, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))]
llps[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
llpsr = data.table(do.call(rbind,lpsr) )
llpsr[, part := factor(part, levels=rev(c("Midpiece","Flagellum")))]
llpsr[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
# motility
lv = list()
lvp =list()
# VAP
m = lm(scale(VAP) ~ Morph, a)
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
lv[['VAP']]=data.frame(response='Average-path (VAP)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lm(VAP ~ Morph, aa)
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$pred <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$motility = 'Average-path (VAP)'
lvp[['vap']] = data.table(newD)
# VSL
m = lm(scale(VSL) ~ Morph, a)
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
lv[['VSL']]=data.frame(response='Straight-line (VSL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lm(VSL ~ Morph, a)
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$pred <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$motility = 'Straight-line (VSL)'
lvp[['VSL']] = data.table(newD)
# VCL
m = lm(scale(VCL) ~ Morph, a)
#summary(m)
#plot(allEffects(m))
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
lv[['VCL']]=data.frame(response='Curvilinear (VCL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
# get predictions
m = lm(VCL ~ Morph, a)
bsim = sim(m, n.sim=nsim)
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
newD=data.frame(Morph = unique(a$Morph)) # values to predict for
X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here
newD$pred <-(X%*%v)
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
predmatrix[predmatrix < 0] <- 0
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
#newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
newD$motility = 'Curvilinear (VCL)'
lvp[['VCL']] = data.table(newD)
llv = data.table(do.call(rbind,lv) )
llv[, response := factor(response, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))]
llv[effect == '(Intercept)', effect:='Independent\n(Intercept)']
llv[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
llv[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
llv[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))]
llvp = data.table(do.call(rbind,lvp) )
llvp[, motility := factor(motility, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))]
llvp[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))] llll = rbind(ll,lls)
llll[, unit := factor(unit, levels=rev(c("single sperm", "male average")))]
llll[effect == '(Intercept)', effect:='Independent\n(Intercept)']
llll[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
llll[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
llll[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))]
g =
ggplot(llll, aes(y = effect, x = estimate, col = response, shape = unit)) +
geom_vline(xintercept = 0, col = "grey30", lty =3)+
geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
#ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.4)) +
#scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
#scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) +
scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) +
scale_shape(guide = guide_legend(reverse = TRUE)) +
#scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
labs(y = NULL ,x = "Standardized effect size",title = 'Sperm-part specific models') +
#ylim(c(0,100))+
#coord_flip()+
theme_bw() +
theme( #legend.position ="right",
plot.title = element_text(size=7),
legend.title=element_text(size=7),
legend.text=element_text(size=6),
##legend.spacing.y = unit(0.1, 'cm'),
legend.key.height= unit(0.5,"line"),
#plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(colour="black", size = 7),
axis.text.y=element_text(colour="black", size = 7),
axis.title=element_text(size=9)
)
g
ggsave(here::here('Output/morpho_effectSizes_virid.png'),g, width = 12, height =10, units = 'cm')
#ggsave('Output/morpho_effectSizes_pair.png',g, width = 12, height =10, units = 'cm')!!! Satellites seem to have the longest tails (and flagellums), faeders the longest midpieces, which translates into the differences in relative measurements. !!!
llcv[effect == '(Intercept)', effect:='Independent\n(Intercept)']
llcv[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
llcv[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
llcv[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))]
g =
ggplot(llcv, aes(y = effect, x = estimate, col = response)) +
geom_vline(xintercept = 0, col = "grey30", lty =3)+
geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
#ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.4)) +
#scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
#scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) +
scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) +
#scale_shape(guide = guide_legend(reverse = TRUE)) +
#scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
labs(y = NULL ,x = "Standardized effect size",title = 'Sperm-part specific models on male CV') +
#ylim(c(0,100))+
#coord_flip()+
theme_bw() +
theme( #legend.position ="right",
plot.title = element_text(size=7),
legend.title=element_text(size=7),
legend.text=element_text(size=6),
##legend.spacing.y = unit(0.1, 'cm'),
legend.key.height= unit(0.5,"line"),
#plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(colour="black", size = 7),
axis.text.y=element_text(colour="black", size = 7),
axis.title=element_text(size=9)
)
g
ggsave(here::here('Output/morpho_effectSizes_virid_CV.png'),g, width = 12, height =10, units = 'cm') llv[effect == '(Intercept)', effect:='Independent\n(Intercept)']
llv[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
llv[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
llv[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))]
g2 =
ggplot(llv, aes(y = effect, x = estimate, col = response)) +
geom_vline(xintercept = 0, col = "grey30", lty =3)+
geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
#ggtitle ("Sim based")+
geom_point(position = position_dodge(width = 0.4)) +
#scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
#scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) +
scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) +
#scale_shape(guide = guide_legend(reverse = TRUE)) +
#scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
labs(y = NULL ,x = "Standardized effect size",title = 'Velocity-type specific models on male velocity') +
#ylim(c(0,100))+
#coord_flip()+
theme_bw() +
theme( #legend.position ="right",
plot.title = element_text(size=7),
legend.title=element_text(size=7),
legend.text=element_text(size=6),
##legend.spacing.y = unit(0.1, 'cm'),
legend.key.height= unit(0.5,"line"),
#plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(colour="black", size = 7),
axis.text.y=element_text(colour="black", size = 7),
axis.title=element_text(size=9)
)
g2
ggsave(here::here('Output/motility_effectSizes_virid.png'),g2, width = 10, height =8, units = 'cm')!!! Against prediction, faeders might have > variable & slowest sperm !!!
Predictions
g1 =
ggplot(a, aes(x = Morph, y = Length_avg)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1.5)+
geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) +
geom_errorbar(data = llp, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
geom_point(data = llp, aes(x = Morph, y =Length_avg), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
#stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_color_viridis(discrete=TRUE)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
labs(title = 'Predictions from sperm part specific models\non male means') +
ylab('Length [µm]') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
#axis.text.x = element_blank(),
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #,
g2 =
ggplot(b, aes(x = Morph, y = Length_µm)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1)+
geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) +
geom_errorbar(data = llps, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
geom_point(data = llps, aes(x = Morph, y =Length_µm), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
#stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_color_viridis(discrete=TRUE)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
labs(title = 'Predictions from sperm part specific models\non single sperm') +
ylab('Length [µm]') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #,
grid.draw(rbind(ggplotGrob(g2), ggplotGrob(g1), size = "last"))
#grid.arrange(g1,g2)
gg1 <- ggplotGrob(g1)
gg2 <- ggplotGrob(g2)
ggsave(here::here('Output/morpho_predictions+boxPlots.png'),rbind(gg2,gg1, size = "last"), width = 15, height =15, units = 'cm') g1 =
ggplot(br, aes(x = Morph, y = Length_rel)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1)+
geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) +
geom_errorbar(data = llpsr, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
geom_point(data = llpsr, aes(x = Morph, y =Length_rel), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
#stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_color_viridis(discrete=TRUE)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
ylab('Relative length to total sperm length') +
labs(title = 'Predictions from sperm part specific models\non single sperm') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #,
g2 =
ggplot(ar, aes(x = Morph, y = Length_rel)) +
geom_dotplot(binaxis = 'y', stackdir = 'center',
position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1.5)+
geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) +
geom_errorbar(data = llpr, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
geom_point(data = llpr, aes(x = Morph, y =Length_rel), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
#stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
scale_color_viridis(discrete=TRUE)+
facet_wrap(~part, scales = 'free_y', nrow = 2)+
labs(title = '\non male means') +
theme_bw() +
guides(x = guide_axis(angle = -45)) +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
plot.title = element_text(size=9)
) #panel.spacing.y = unit(0, "mm")) #,
gg1 <- ggplotGrob(g1)
gg2 <- ggplotGrob(g2)
grid.draw(cbind(gg1, gg2, size = "last"))
#grid.arrange(g1,g2)
ggsave(here::here('Output/morpho_predictions+boxPlots_rel.png'),rbind(gg1,gg2, size = "last"), width = 7.5, height =15/4, units = 'cm')
# End